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Abstract. This paper deals with a high-order accurate implicit finite-difference 
approach to the pricing of barrier options. In this way various types of barrier 
options are priced, including barrier options paying rebates, and options on 
dividend-paying-stocks. Moreover, the barriers may be monitored either con- 
tinuously or discretely. In addition to the high-order accuracy of the scheme, 
and the stretching effect of the coordinate transformation, the main feature of 
this approach lies on a probability-based optimal determination of boundary 
conditions. This leads to much faster and accurate results when compared 
with similar pricing approaches. The strength of the present scheme is partic- 
ularly demonstrated in the valuation of discretely monitored barrier options 
where it yields values closest to those obtained from the only semi-analytical 
valuation method available. 



1. Introduction 

Barrier options are a type of path-dependent options whose values depend on 
the specific path followed by the underlying asset during the option's life, and w.r.t. 
some specified asset values, usually referred to as barriers. These options are exotic 
derivatives traded in over-the-counter markets and they can therefore be tailored to 
specific customer needs. Moreover, the additional constraints imposed on them by 
the barrier makes them a much cheaper and attractive product then the standard 
options. Barriers can also be added to any existing type of standard or exotic 
option. The barrier option market has thus been expanding rapidly and it has 
been estimated that it has doubled every year since 1992. In parallel with this 
development, a large number of newly designed barrier options have been trading 
very actively in financial markets [71 130j . 

Analytical formulas exist for most of the standard barrier options. Rubinstein 
and Reiner [28], Rich [26], and Kunitomo and Ikeda [21] derived analytical formu- 
las for a variety of standard knock-in and knock-out European options with full 
barriers. Heynen and Kat jl8j obtained similar formulas for some special types of 
barrier options, namely the partial barrier options, and for so-called outside bar- 
rier options where it is another variable different from the underlying asset which 
determines whether the option knocks in or out. Closed- form solutions for double 
barrier options, some of which are time-dependent, have been obtained [2T | [TS t [25] . 
The extension of these formulas to the more complex case of American-style bar- 
rier options has been undertaken by Broadie and Detemple |6j, Gao et al. [14], as 
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well as Haug ^16]. However, all the formulas obtained in such extensions are only 
closed- form approximations. 

Despite all these efforts to derive analytical pricing formulas for barrier options, 
there is still no analytical formulas available for large classes of both European- 
and America-style barrier options [16l [21]. Although Merton |23j proposed the 
closed-form solution for the continuously monitored down-and-out barier option 
in 1973, it's only in the 1990's that valuation formulas were obtained for other 
variants of this European-style option [551 [Ml [HI HI] • This is simply because the 
valuation of financial options has led to mathematical models which are most often 
challenging to solve. For instance, for most of the complex options, there are no 
analytical formulas available, and almost all formulas have been obtained under 
the assumptions that, amongst others, the underlying asset price has a lognormal 
distribution with constant drift and volatility parameters, and that there are no 
transaction costs [551 [TB]. However, as empirical evidence suggest the contrary, new 
models with stochastic or time-dependent volatilities, or including transaction costs 
have been designed and implemented. Nevertheless, these attempts to improve on 
the assumptions underlying the derivation of analytical formulas have been made 
mostly for standard American options [TTI [HI [20] . 

For the valuation of options, and exotic options in particular, numerical methods 
remain a tool of choice. This is first justified by the fact that the number of exotic 
options that can be designed is limitless, and as new products enter the market one 
of the most practical way, if not the only available means to price them, are quite 
often the numerical methods. They are also a convenient benchmark for testing 
the validity of analytical formulas, especially as they become available. There are 
even instances where they are faster than analytical methods, even those analytical 
methods involving fast converging infinite sums, when the required accuracy isn't 
too high. 

As a numerical method for option valuation, the lattice methods have been used 
extensively, and for barrier options it appears however that they are essentially 
useless without the proper positioning of the barrier which must line-up with the 
tree nodes. Despite all the techniques that have been developed to improve on 
lattice methods, they are still quite computer- intensive ,16J. Monte Carlo Methods 
have been considered as inflexible and unreliable for option pricing, until the last 
decade where they yielded more promising results based on innovative techniques 
which are now a topic of current research. Barraquand and Martineau [5] amongst 
others, obtained in this context some interesting results for multi-asset American 
options. The most commonly used numerical method in option pricing these days 
appears to be the finite difference methods. They have an acceptable computational 
cost when an appropriate implicit method is used (see [301 [22l [20] ) . They can also 
be extended to cater for American-style options, by formulating them as linear 
complementary problems. 

In this paper we consider an implicit finite difference approach to the valuation of 
barrier options. These barrier options may have various features including double 
barriers, dividend-paying-stocks, rebate payments, or some combination of these. 
Barriers may be monitored either continuously or discretely. We use a well-known 
0-metliod which is fourth-order accurate in space and second-order accurate in 
time, by relating the parameter 6 to the mesh sizes. In addition to the high-order 
accurate scheme and the stretching effect of the coordinate transformation used. 
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the strength of this approach hes on a probabihty-based optimal determination of 
the boundary conditions, along which the option values are known exactly. In this 
way, with a reasonable accuracy and for the most practical considerations, we are 
generally able to use fewer than 20 asset prices and 20 time steps per year to value 
most of the barrier options, especially when the barrier is continuously applied. 
However, the efficiency of our approach is much clearly demonstrated in the case 
of discretely monitored barrier options, when we compare the results we obtained 
with that of several other authors. As usual, due to some parity considerations, we 
shall focus our attention only on knock-out call options. 

This paper is organized as follows. In Section [2] we discuss the mathematical 
model for valuing barrier options. Section [3] is devoted to the discretization of the 
modeling differential equation, while Section[4]discusses applications of the resulting 
implicit scheme to the various types of barriers options, and Section [5] reports the 
corresponding numerical results. 



2. Option pricing model 

We make the usual assumption that the underlying asset price S has a geometric 
Brownian motion, with drift and volatility parameters fi and cr, respectively. Thus 
the process followed by S can be written with the usual notation and in terms of 
the Wiener process W as 

dS^fidt + adW. (2.1) 

The price / — f{S,t) of a European option contingent on S must then satisfy 
the Black-Scholes differential equation 

C-f^O (2.2) 
where the partial differential operator £ is given by 

and where r denotes the risk-free interest rate. We shall usually assume that the 
underlying asset price pays a continuous dividend yield at a continuous rate of q 
per year, and as we are placing ourselves in the traditional risk-neutral world in 
which the market price of risk is zero, we shall have in this case fi — r ~ q. 

Thanks to the put-call parity for European options, and also to the similar parity 
relation between knock-ins and knock-outs, wc shall consider only knock-out call 
options. To begin with, let B denote the constant barrier value for a down-and- 
out barrier option. Since the value of the option at expiration time T is its payoff 
Max(S' — K,Q), where K is the strike price of the option, the initial condition for 
equation (12. 2|) is given by 

/(S',T) =Max(S'-X,0), forS'>S. (2.4) 

The corresponding boundary conditions follow from the properties of call options. 
In the simplest case where the constant barrier value B is continuously applied, 
denoting by Rb the rebate received if the barrier is ever breached, the boundary 
conditions are 
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f{B,t)^Rb, yt (2.5) 
f{S,t)r^S, yt, as 5*^00. (2.6) 

It should be noted that boundary condition (|2.6p which is usuaUy used for the 
valuation of call options is not suitable for an optimal algorithm, in particular 
because the value of S in that condition is generally exceedingly large, and the 
resulting equality is only an approximation, and all these tend to slow down the 
valuation algorithm. We shall therefore always strive in this paper to find optimal 
boundaries of the solution domain along which the option values are known exactly, 
and which also yield an optimal truncation of the solution domain. 

The Black-Scholes equation (|2.2p can be analyzed directly for the valuation of 
barrier options, as for example in [4] and |30j . However, there is a numerically more 
effective coordinate transformation that we shall adopt here, and which is obtained 
by first setting 

X = ln{S/K), (2.7) 

r=y(T-t). (2.8) 

Since for all practical considerations the underlying price range generally lies in a 
subinterval of {K/3,3K) , it appears that (|2.7p restricts the underlying asset price 
in the a:-coordinate to range most often in a small interval such as (—1.099, 1.099) . 
On the other hand, (|2.8p not only transforms (|2.2p into a forward-time problem, 
but it also has the advantage of shrinking the time range. Under the change of 
coordinates (|2.7p and (|2.8p . equation (|2.2p is transformed into 



— = — + (^^-1)— -J^i./ 
dr dx"^ dx 

where v = 2/i/(T^, and vi = 2r/a'^. So, v = vi — V2, where V2 = ^qja^ . Here v\ and 
vi are the only two non-dimensional parameters of the Black-Scholcs equation for 
dividend-paying-assets. To discard the terms in / and df /dx in this last differential 
equation, we now make a change of the dependent variable by setting 

f{S,t)^Ke'^-+'<^u{x,T), (2.9) 

where a = —j{v—\) and 7 = — i(i' + l)^ — 1^2- This last change of variables reduces 
the original equation (|2.2p to the diffusion equation 

du d^u 

d-r = dx^- _ ^ ^'-''^ 

Letting Xh = \tl{B / K)^ the initial and boundary conditions given earlier are also 
transformed into the following expressions in the [x, r)-coordinates. 

m(x,0) = Max(e(i/2)(i^+i). _g(i/2)(.-i)a,^0^ ^ ^ > (2.11a) 
u{x, t) ~ e^i-")^-''^, s&x^oo (2.11b) 



and 



u(a;b,r) = e-(""''+^"). (2.11c) 
K 
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It is worth nothing at this point that the change of variable (|2.9p is also very 
convenient for an efficient numerical algorithm. Indeed, it gives the final solution 
as a scalar multiple of the numerically computed factor u(x,t), where the scalar 
j^^ax+jT known exactly and so, does not involve any error. In this way any error 
expansion due to transformation (j2.9[) will tend to be minimized. 



3. Finite difference scheme 

Denote hy h — Ax and k — At the mesh sizes in the space and the time 
directions, respectively. Denote also by C/" the discrete approximation of u{x,t) 
at the grid node {xj,tn), where xj = jAx, and i„ = nAt. Let S'^ be the difference 
operator given by <5^J7" = Uj^-i — 2?7" + f/j+i- For the discretization of the reduced 
equation (|2.10p we choose the two-time level, three-space-point scheme 

jjn+l _ jjn , 

' A, ' = {(^5lur' + (1 - ^) ^f) (3-1) 

also called weighted average approximation or (^-method. The weight 6 in this 
expression is assumed to satisfy the condition < < 1 , in order to avoid negative 
weights. The values 6 = and 9 = 1 give the explicit scheme and the fully implicit 
scheme, respectively. If we let /3 = At/(Aa;)^, i.e. /3 = k/h? , then for < < |, 
the scheme is stable if and only if /3 < 1/2(1-26'). For \ < < 1, it is 

stable for all values of /3. For this scheme, the truncation error, i.e. the amount by 
which the exact solution of the differential equation does not satisfy the discrete 
approximation (13. 1|) . is generally of first order accuracy in At, while the popular 
Crank-Nicolson scheme which corresponds to 9 = 1/2 is second-order accurate in 
both At and Aa;. 

For a much higher accuracy of the finite difference scheme (fds) (|3.ip . we shall 
relate in this paper the choice of the parameter 9 to the values of Ax and At of the 
grid line spacings. More precisely, we shall assume that 

It is well-known that the resulting scheme is O ((Ai)^ -t- (Aa;)*) , i.e. second- order 
accurate in time and fourth order accurate in space (see e.g. [H]). Such a scheme 
has been used in [22] for the valuation of American options with error correction at 
the critical boundary, leading to results that are comparatively more accurate than 
the standard ones found in the literature. However, this specific weighted average 
scheme has not yet been applied to the valuation of barrier options, to the best 
of our knowledge. Zvan et al used a method called point-distributed finite volume 
scheme in [30] to directly discretize the Black-Scholes equation, but that method is 
only first order accurate in time, and the paper does not comment on the speed (in 
terms of number of time steps) of the scheme for achieving a given level of accuracy. 

Since /3 is positive, condition p.2p restricts 9 to the interval [0,1/2), but the 
condition for stability /3 < 1/(2(1 — 29) is always satisfied in this case. Moreover, 
the condition < 9 implies that 

(Ax)^ < 6At. (3.3) 

This in practical terms means that even by requesting (|3.2p to hold, we can still 
take large time steps while maintaining accuracy and stability. 
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If we denote by M and L the number of space- and time-intervals, respectively, 
in the solution grid, and by Wh.k the finite difference operator given by 

Wh.k ■ = -mf-i + (1 + wo)uf - m^+i. 

then an expansion of the discretized equation p.ip leads to a linear algebraic system 
of equations of the form 



WH.k-Uf^^ = Wh.k-Uf, forn = 0,...,L-l andj-l,...,M-l. (3.4) 

The system p.4p can be written in terms of tridiagonal matrices and we shall 
solve it with the Thomas algorithm which includes Gaussian elimination without 
pivoting. 

4. Numerical solution for barrier options 

We discuss in this section the application of the fds described above to the 
derivation of a numerical solution to the pricing problem for barrier options. 

4.1. Overview. As already indicated, we shall always implement whenever pos- 
sible in this paper an optimal truncation of the solution domain by replacing the 
approximate boundary condition of the form ()2.6p . or (j2.11b|) equivalently, with 
a boundary condition for which both the option value is known exactly and the 
resulting solution domain is the smallest possible. Let 5*0 be the initial price of the 
underlying and S = St the price of the underlying at a time t, < i < T. In the 
case of a down-and-out call option, the probability P{St > B) that the barrier will 
not be breached at time t is given by 

P{St >B) = $(a*), (4.1) 

where $ is the standard normal distribution function and 

a* = (ln(5'o/B) + fi*t)/aVi, with ^i* = fi - 

Let 6 be the smallest number such that ^{S) has numerical value 1. The value 
of S depends on the level of accuracy required and also on the computing system 
used for the calculations. An appropriate value for 6 usually lies in the interval 
(3.7,6.5), depending on the level of acuracy required. For a given parameter set, 
the minimum value of 5*0 which guarantees that the barrier will not be breached 
at time t is given by the inequality a* > S, or equivalently by Sq > Be^'^^~^^ *. 
Now, let tp be the turning point of the function t i—s- 8a\ft — ij,*t, and denote by Sm 
the minimum value of the current stock price So above which the barrier becomes 
worthless, and set Xm = ln(5',„/if). Then it is easy to see that 

_ j xb + SaVT - ^i*T, if /i* < or > T 
I Xfc + ioyftp — /i*ip, otherwise. 

For any value of Sq above Sm, or equivalently for any value of xq above Xm, 
where xo = ln(5o/-ftr), the barrier is worthless and the option behaves exactly as a 
standard option. In particular, the corresponding option value along the boundary 
S — Sm or above this boundary is given by the well-known Black-Scholes formula 
for European call options. 

It therefore follows from Equations (|2.5p . (12. 8p and (|4.2p that the solution domain 
for u{x, r) reduces to the rectangle [xb, Xm] x [0, ^T]. Similarly for an up-and-out 
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barrier option, if we let Sm denote the maximum value of the current underlying 
asset price below which the barrier becomes worthless, then Xm — ln(<S'm/if) is 
given by 

j Xb ^ SaVT - fi*T, iifi* >Ooitp>T ,^ , 

x„i = < (4.3) 
I Xh — 8a^/T^ — ^*tp, otherwise. 

In this case the optimal solution domain reduces to the rectangle [a;T„,a;;,] x 
[0, ^T]. This optimal determination of the boundary conditions has several appli- 
cations, including the classification of double barrier options. A similar determina- 
tion of optimal boundary conditions can be achieved for time- varying barriers, as 
opposed to the constant barriers we have discussed here. Note that by construction 
we have Xm > Xb in (|4.2p and a;,„ < Xb in (|4.3p 

Our determination of the numerical solution for a pricing problem will amount 
to calculating u{xq, tm), where tm = and this may require interpolation if xq 
does not line up with a grid line. However, for M large, it is always possible to 
slightly enlarge the range of the spacial variable x without affecting the performance 
of the algorithm in any way, just by adjusting the value of Xm by a fraction of the 
grid increment By so doing, we can always line up xq with a grid point, and such 
adjustment of Xm is feasible when for instance the grid increment h = Ax is small. 
We shall therefore combine these two techniques in our calculation of u{xo,tm)- 

The price of barrier options often displays some singularities near the barrier 
value and it is therefore customary in a finite difference pricing approach for these 
options to use adaptive grids in an attempt to achieve the required accuracy of 
the numerical solution with a smaller number of mesh points, i.e. with a faster 
algorithm. Suppose that the grid increment h ~ Ax of the spacial variable x is 
non constant and denote by hj — Sj+i — Xj the jth grid increment, where the x^s 
are grid points in the direction. The number qj — hj/hj-i is often called non 
uniformity parameter of the difference grid [T]. A direct substitution of the non 
constant increment hj into the FDS (jS.ip gives rise to the nonuniform version of 
equation (|3.4p . A similar but much simpler discretization of equation (j2.10p can be 
derived for a nonuniform time stepping. For the nonuniform discretization of either 
of the variables, we used a geometric progression. Thus for the space variable, we 
set hi = ft,, hj = W^^h for some positive numbers R and h. Here, R turns out to 
be the non uniformity parameter and the ft^'s form a geometric progression with 
common ratio R. Values of R not close to 1 will lead to instability of the numerical 
solution, while for R close to 1 the resulting scheme will tend to be stable [Tj . 

A two-dimensional grid made up of M space intervals and L time intervals will 
be referred to as a grid of mesh (size) M x L. In the case of the fds that we are 
considering, L and M are to be so chosen that (|3.3p is satisfied, and this is easily 
achieved by letting L = M. In some very rare cases (e.g when tr > 35%), the 
occurrence of negative option prices at the intermediary time steps might become 
persistent and slow down the convergence of the scheme. This type of difficulty is 
also easily resolved by letting L be slightly larger than M, say L = 1.5M. 

It should be noted that our algorithm for the implementation of the optimal 
boundary condition shares some similarities with the method for finding the con- 
ditional expectation estimator in the Monte Carlo simulation of a down-and-in 
barrier option [27j . Indeed, in such a Monte Carlo valuation, as soon as the barrier 
is breached, the simulation run is ended, and the corresponding estimator is the 
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Black-Scholes value for the resulting option parameters. In the same way for the 
finite difference scheme, as soon as we are certain that the initial price is sufficiently 
far away for the barrier to become worthless, the corresponding option value is its 
Black-Scholes value. 

4.2. Discretely monitored barriers. In the case of discretely sampled barriers, 
we assume that barrier monitoring occurs either on a daily basis or on a weekly 
basis. We adopt a day count convention which will allow a comparison of our nu- 
merical results with others at our disposal by letting a year consist of 50 weeks, 
and a week to consist of 5 days, so a year has 250 days. For the numerical solution, 
we let p be the number of time steps implemented between two consecutive moni- 
toring dates, so that the total number of time steps is given by L = iVp, where N 
is the number of monitoring dates during the option's life. As is well-known [H [8] , 
for increased accuracy and improved convergence, we shall place the barrier itself 
midway between grid points. 

5. Numerical results 

We start by showing the result of the performance test for the proposed optimal 
truncation of the solution domain. In most papers where approximate conditions of 
the form (j2.6p are used to determine boundary conditions in barrier option pricing 
problems, the value used for Smax, the maximum asset price in the solution domain, 
is usually not indicated. That is the case for instance in [30] and [4]. However, by 
rewriting and implementing the same algorithm described in the latter paper under 
the scheme termed modified explicit finite difference (mefd) in that paper, we've 
found by trial and error that choosing S'max = 2S'o -I- 200 gives exactly the same 
results as those in Table 1 and Table 4 of the said paper. This value of S'max also 
turns out to give the best results for the parameter sets considered in that same 
paper. Other choices for S'max which are commonly used for a similar range of data 
are Smax = 2So, or Smax — Sq -\- 100, as in [5] and [52]. However, inappropriate 
choices for Smax lead to very poor results as we show in Table [1] 

Contrary to conditions of the form given by equations (|4.2|) - (l4.3p that yield a 
systematic determination of the optimal boundary conditions, there is certainly no 
formula available for the best choice of Smax and it's value is usually determined 
only by trial and error. For two different values of Sq, Table [T] compares for a 
down-and-out barrier option the convergence rates under the mefd scheme and 
under our optimal boundary explicit scheme (obes ) using data from Table 1 of [4] . 
The two different formulas used for Smax are Smax = 2So-l-200 and Smax = 2So. The 
parameter set is T = l,K = 100, B = 90, ct = 0.25, and r = 0.10. For the explicit 
scheme used, the mesh size is determined by a single parameter L representing the 
number of required time steps, because of the usual stability conditions relating 
the time and the space increments in an explicit scheme. Moreover, this constraint 
for the MEFD scheme upon which all explicit schemes considered in this paper are 
based is given by 

Ay = ACT\/At, (5.1) 

where y — log(S). Numbers within parentheses in Table [1] are CPU times. All runs 
in this paper were carried out on a conventional Pentium Celeron 2.40 Ghz CPU. 
In all the tables the parameters q and Kb are assumed to be zero, unless otherwise 
indicated. 
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Table 1 . Comparison of the accuracy and the CPU time under the optimal 
boundary explicit scheme (OBEs) and the MEFD of [4] for a down-and-out call. 
The fixed parameter set is T = I, K = 100, B = 90, = 0.25, r = 0.10. 
Numbers within parentheses are CPU times. 



Mesh So = 95 

(L) Closed-Form: 5.9968 



So = 91 

Closed-Form: 1.2738 



OBES 



MEFD 



OBES 



MEFD 



Sma.K — 

25*0 + 200 



2^0 



25o + 200 



50 


6.01109 


6.0111 


6.8299 


1.2654 


1.2654 




(0.36) 


(0.344) 




(0.39) 


(0.33) 


100 


5.9984 


5.9984 


6.8603 


1.2704 


1.2704 




(0.53) 


(0.44) 




(0.53) 


(0.47) 


150 


5.9965 


5.9997 


6.7449 


1.2719 


1.2719 




(2.09) 


(2.33) 




(2.203) 


(2.42) 


200 


5.9993 


5.9993 


6.6897 


1.2725 


1.2726 




(2.54) 


(2.70) 




(2.58) 


(2.85) 


700 


5.9970 


5.9970 


6.74 


1.2738 


1.2738 




(15.81) 


(17.766) 




(15.75) 


(18.16) 


800 


5.9972 


5.9972 


6.78 


1.2739 


1.2739 




(19.20) 


(21.906) 




(19.44) 


(22.30) 


900 


5.9973 


5.9973 


6.80 


1.2739 


1.2739 




(10.703) 


(12.047) 




(10.922) 


(11.625) 


1000 


5.9969 


5.9969 


6.7991 


1.2739 


1.2739 




(33.86) 


(40.59) 




(34.31) 


(42.25) 


2000 


5.9970 


5.9970 


6.7949 


1.2738 


1.2738 




(72.67) 


(85.97) 




(75.12) 


(85.797) 


3000 


5.9969 


5.9969 


6.7718 


1.2738 


1.2738 




(189.34) 


(238.141) 




(192.062) 


(236.66) 


4000 


5.9969 


5.9969 


6.8020 


1.2738 


1.2739 




(263.28) 


(320.187) 




(268.39) 


(335.31) 



Columns 2 and 3 of Table [T] show that the OBES and the mefd schemes give es- 
sentially the same option values for the given parameter set, when 5niax = 250 + 200 
in the mefd. However, our OBES appears to be up to 20% faster for all significant 
number of time steps (i.e. for L > 150). For small values of L below 150, the excess 
amount of time used in the OBES to calculate the option values at different time 
steps on the upper boundary determined by Sm exceeds the excess time required in 
the mefd to implement the additional time steps. But this is of no importance as 
regards the performance of the schemes, as for explicit schemes accuracy is generally 
not achieved below L — 1000 steps. The same conclusion applies when comparing 
columns 5 and 6 of the table. 
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A comparison of Column 3 and 4 of the same table shows the importance of 
choosing an appropriate value for S'max- Indeed, the same calculations in Column 
3 when done in Column 4 with S'max = 25*0, which as already mentioned is a 
commonly used value for similar parameter sets in option pricing problems, reveals 
that the scheme certainly does not even converge to the closed- form value of 5.9968. 
However, even for those values of S'max for which the mefd does converge, there is 
no systematic way for finding them. 

Table 2. Effect of volatility on the performance of the OBES 
and the MEFD of |4] for a down-and-out call with fixed parameters 
50 = 95, T = 1,K = 100,B = 90, and r = 0.10. 



Closed- 
Mesh OBES MEFD form 



•S'max "S'max 

= 25o + 200 = 2So 



a = 


0.25 


200 


X 


200 


5.9993 


5.9993 


6.6897 


5.9968 






2000 


X 


2000 


5.9970 


5.9970 


6.7949 




a = 


0.30 


200 


X 


200 


5.9071 


5.9075 


7.6180 


5.9060 






2000 


X 


2000 


5.9061 


5.9065 


7.7124 




a — 


0.40 


200 


X 


200 


5.7498 


5.7802 


9.2320 


5.7502 






2000 


X 


2000 


5.7503 


5.7789 


9.5294 





The validity of our optimal determination of boundary conditions is further 
demonstrated by the effect of volatility in Table [21 Indeed, this experiment shows 
that as the volatility increases, the MEFD becomes more and more inaccurate for 
almost every possible choice of S'max, and does not even tend to converge for values 
of a above 40%. On the other hand the OBES is perfectly unaffected by any change 
in the value of a. Indeed, for large values of cr, the approximation f(S,t) ~ S for 
all t is only tolerated in the MEFD for much larger values of S. We've found that 
by adequately increasing the value of Smax for a larger value of cr, convergence 
under the MEFD can still be achieved, but at the expense of much longer computing 
times and additional trial and error experiments to find S'max- Some complications 
might also arise in the choice of S'max, due to the constraint ()5.1|) . As for the OBES, 
any change in the value of a is fully captured by both S„i and the corresponding 
Black-Scholes value, so that the OBES remains totally unaffected by such changes. 

The efficacy of combining the high-order accurate implicit scheme determined by 
equation (|3.2p . with the stretching coordinates transformation given by (|2.7p - (|2.8p . 
and the optimal truncation of the solution domain is first tested for a down-and-out 
barrier call option in Tabled We denote by hobis (high-order optimal boundary 
implicit scheme) the high-order implicit finite difference scheme that applies the 
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Table 3. Down-and-out option values for initial asset prices closest to the 
upper boundary (U), the lower boundary (L) or at midway (M). The fixed 
parameters are K = 150, and B = 180. Numbers within parentheses are 
CPU times. 





T = 0.25, (T = 0.20, r = 1, (T = 0.45, 
r = 0.05 r = 0.07 


= 271.906 Srri = 1345.08 
Xm — Xb — 0.4125 x„i — Xb — 2.0112 


So Closed- HOBis HABis So Closed- hobis 
Form Mesh Mesh Form Mesh 


u 


271.905 123. 7d8 3x1 25 x 25 1345.07 1205.21 2 x1 
(0.00) (0.05) 

271.902 123.765 3 x 1 25 x 25 1345.00 1205.14 3 x 1 
270.000 121.862 10 x 10 25 x 25 1344.00 1204.14 7 x7 
265.000 116.861 12 x 12 25 x 25 1340.00 1200.14 10 x 10 


L 


180.001 0.0026 6 x 6 30 x 30 180.001 0.0015 3 x 1 
(0.03) 

180.010 0.0264 11 X 11 30 x 30 180.01 0.0015 3 x 1 
181.000 2.6297 135 x 135 400 x 400 181.00 1.4817 95x 95 
(1.20) (14.25) (0.55) 


M 


225.953 77.2335 95 x 95 600x600 762.54 622.632 85 x 85 
(0.53) (39.26) (0.52) 



optimal determination of the boundary conditions, while the corresponding scheme 
in which a boundary condition is determined by the approximate condition p.6p 
will be referred to as habis (high-order approximate boundary implicit scheme). 

A striking feature of the high-order accurate hobis is its fastest convergence for 
initial asset prices closest to boundaries of the solution domain. For initial asset 
values closest to 5*^, on the upper boundary (U), the minimum mesh size that 
returns the anlytical option price is 2 x 1 and corresponds to a CPU time of 0.00 
seconds. For values closest to the barrier on the lower boundary (L), the minimum 
mesh size is 7 x 7 and corresponds to a CPU time of 0.02 seconds. This performance 
is clearly due to the exactness of option values at the boundaries in combination with 
the fact that interpolation yields best results for points closest to known function 
values. But this is also a consequence of the coordinates transformation (in the 
first set of data with a = 0.20 in Tabled Xm — Xb = 0.4125 < 1), the higher-order 
accuracy of the scheme that forces all error less than one to quickly disappear, and 
the implicit nature of the scheme that also tends to quickly reduce any error. For 
instance, in the second set of data in Table [3] with a = 0.45, and So = 1345.00, 
although the space increment is 



Ax = 2.011/2 = 1.005 > 1, 
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all significant errors have disappeared by the completion of the first time step. The 
table also shows as expected that the more we move away from the boundaries, 
the higher the error tends to become in option values with a fixed mesh size. We 
have therefore also indicated the minimum CPU time required for initial values 
midway (M) between the boundaries, and we can reasonably argue that for any 
given scheme and any parameter set, the maximum CPU time for all values of Sq 
in the solution domain is close to the maximum CPU time obtained for Sq close 
to the boundaries and at midway. This means for instance that for the first set 
of data with a — 0.20, and the second set of data with a = 0.45, to within four 
significant digits of accuracy, the maximum CPU time is close to 1.20s and and 
0.55s, respectively. It should be noted that this is also the maximum CPU time 
for any possible value of the initial price. Indeed, another feature of the optimal 
boundary condition is to determine exactly when a barrier becomes worthless for a 
given value of 5*0, so its value is reduced to the corresponding vanilla Black-Scholes 
value. More precisely, by construction of Sm, whenever Sq > 5*^ the barrier option 
value is the vanilla Black-Scholes value, and so any scheme (explicit, implicit or 
otherwise) based on the optimal boundary returns the option value in zero seconds 
for such values of Sq. 



Figure 1. Convergence patterns under the HO- 

BIS for a double knock-out. The parameter set is 

So = 100, T = 0.5, K = 100, Bi = 75, Bu = 125, cr = 0.20 
and r = 0.10. All absolute errors are less than 10% with a 17 X 17 mosh, and 
less than 2 X 10"'' with a 70 X 70 mesh. 




The 5th column of Table[3]shows that implementing the same high-order accurate 
implicit scheme by using the approximate boundary condition of the form (|2.6p leads 
to much more inaccurate results and requires about 10 times more computing time, 
although the value of 5'max = 250 + 200 appears to be the best choice for all sets of 
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parameters considered. This numerical experimentation confirms the fact that the 
superior accuracy of the hobis is valid for all finite difference schemes. 

The HOBIS curves depicted in Figure [T] for a double knock-out also display a sim- 
ilar trend of convergence described in Table [3l in which accuracy is much quickly 
achieved for points closer to the boundaries of the solution domain, while the middle 
point returns the most coarse value (at least for the first few mesh sizes). The pa- 
rameter set for this figure is 5*0 = 100, T = 0.5, if = 100, = 75, Bu = 125, cr = 0, 
and r = 0.10. The figure clearly shows the global accuracy of the hobis for different 
mesh sizes. For a 17 x 17 mesh size, the hobis curve is essentially indistinguishable 
from the analytical curve and the maximum absolute error is 5.1 x 10~^. This error 
reduces to 3.0 x 10^"^ for a 70 x 70 mesh size. 

Table 4. Maximum absolute error under the HOBIS for continuously moni- 
tored down-and-out and double knock-out calls. The fixed parameter sets are 
K = 100,5 = 90, r = 0.10 for the down-and-out call and K = 100, Bl = 
75, Bu = 125, r = 0.10 for the double knock-out call. The maximum absolute 
error is that on the M + 1 computed option prices for any given mesh size of 
the form M X L. 



DOWN-AND-OUT CALL DOUBLE KNOCK-OUT CALL 





para- 
meters 


mesh 


error 




para- 
meters 


mesh 


error 






T = 0.5 


7x7 


2.4 X 10" 


-2 


T = 0.5 


20 X 20 


7.9 X 10- 


-2 




T = 0.5 


20 X 20 


1.3 X 10" 


3 


T = 0.5 


100 X 100 


3.1 X 10" 


-3 


0.15 


T = 0.5 


100 X 100 


3.0 X 10" 


-5 


T = 0.5 


200 X 200 


7.8 X 10- 


-4 




T = 1 


100 X 100 


3.6 X 10" 


-5 


T = 1 


200 X 200 


4.0 X 10- 


-4 




{T^l, 
Rb = 3} 


100 X 100 


5.6 X 10" 


-5 


T = 0.1 


200 X 200 


3.5 X 10- 


-3 


0.35 


T = 1 


100 X 100 


2.6 X 10" 


-4 


T = 1 


200 X 200 


6.3 X 10- 


-2 



For almost all typical model parameters, the hobis and analytical curves are 
essentially indistinguishable in the case of a down-and-out call, even for the smallest 
3x1 mesh size. We've therefore given in Table |4] the maximum absolute error on 
computed option prices for some given mesh sizes, and for both a down-and-out 
call and a double knock-out call. The table shows that the absolute hobis error 
is strictly less than 10~^ for a 20 x 20 mesh when the option is a down-and-out, 
and that global accuracy is much quickly achieved with a down-and-out than with a 
double knock-out. On the other hand there does not seem to be a simple correlation 
between accuracy and expiry date for a fixed mesh size in the case of the double 
knock-out call, at least for the three different expiry dates considered. However, 
it clearly transpire from the table that the hobis is very sensitive to volatility, 
especially in the case of a double knock-out call where the maximum absolute error 
is about 100 times larger when the volatility goes from 0.15 to 0.35. The table also 
shows that accuracy is a bit more expensive when the computation of option values 
incorporates a rebate payment. 
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Figure 2. Error in the computed option prices / under the HOBIS, the 
Crank-Nicolson, and the fully implicit schemes. The parameter set is T = 
0.5, K = 100, _B = 90, cr = 0.20 and r = 0.10. The mesh size is 40 X 40. 
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Although the hobis is fourth-order accurate in the space variable and second- 
order accurate in the time variable while the popular Crank-Nicolson scheme is 
second order accurate in both the space and the time variables, this does not 
necessarily precludes the actual numerical implementation of these schemes to prove 
otherwise. We've therefore plotted against the initial price 5*0 in Figure [D the 
numerical errors in the computed down-and-out option prices under the hobis, 
the Crank-Nicolson, and the Fully implicit schemes. The latter scheme which was 
used for option pricing in [30] is only first-order accurate in time. The figure 
shows that the hobis is much more accurate than the Crank-Nicolson scheme for 
the corresponding mesh used, and the fully implicit scheme performs very poorly 
compared to the hobis. In fact the maximum absolute errors resulting from these 
plots are 0.00193 for the hobis, 0.00466 for the Crank-Nicolson, and 0.02392 for 
the fully implicit scheme. The error functions for the last two schemes have each a 
singularity where the error tends to jump to zero, but this has essentially no effect 
on their performance. We have however, plotted the same error functions in Figure 
[3] for the two best performing schemes, namely the HOBIS and the Crank-Nicolson 
scheme, to clarify the difference between their levels of accuracy. Although the only 
mesh size used for these two figures is 40 x 40, our experimentations shows that this 
trend persists for all mesh sizes larger than say, 20 x 20. The common parameter 
set for the two figures is T = 0.5, K = 100, B = 90, cr = 0.20 and r = 0.10. 

Table [5] compares hobis values with those reported in Table 1 and Table 2 of 
[50] for a down-and-out and a double knock-out call, and corresponding to values 
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Figure 3. Error in the computed option prices / under the HOBIS, and the 
Crank-Nicolson schemes. The parameter set is T = 0.5, K = 100, B = 90, a = 
0.20 and r = 0.10. The mesh size is 40 X 40. . 
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obtained in [30] using a fully implicit method of first order of accuracy (zvan), 
and in [8] using a trinomial tree method (c-v). The options parameters used are 
5- = 100, T = 0.5, K = 100, a = 0.20, r = 0.10. For the down-and-out, B = 99.9, 
while for the double knock-out, the lower barrier is set at Bi = 95 and the upper 
barrier is set at Bu = 125. For continuously monitored barriers, the hobis values 
correspond to the analytical values to within the required accuracy, as are those 
computed by Cheuk and Vorst, while surprisingly the ZVAN values of [30] do not 
coincide with analytical values. Closed-form solutions for continuously monitored 
barrier options have been derived by various authors [131 UHl (HI [HI US]- The 
hobis appears not only to be more acurate, but also about twice faster than the 
implicit scheme of [SQl- In the case of discrete monitoring, the HOBIS values are 
essentially the same as the C-v values to within the required accuracy, but differ 
significantly again from those obtained using the fully implicit method of [30]. In 
general accuracy is very costly in terms of CPU time for discretely monitored barrier 
options and the computing time in Table [S] for discrete monitoring corresponding 
to the ZVAN column looks too little. For similar amounts of time the HOBIS returns 
similar option values, but does not converge to those values, and it generally requires 
up to a hundred of seconds to converge significantly. 

Although Table [5] clearly demonstrates the higher performance of the hobis in 
terms of convergence and accuracy for continuously monitored options, it is hard 
to make the same conclusion from this table about discretely applied barriers, and 
this is largely due to the non availability of exact solutions in such cases. 

We have therefore compared in Table [6] the hobis values for a discretely moni- 
tored down-and-out option with those taken from Table 2 of [12]. Fusai et al [T2] 
have reduced the valuation problem for discretely monitored down-and-out barrier 
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Table 5. Down-and-out and double knock-out call values with continuous 
and discretely applied constant barriers. 



Monitoring Down-and-out Call Double Knock-out Call 
irequency 



HOBIS ZVAN C-V HOBIS ZVAN C-V 



Continuous 0.165 0.164 0.165 2.033 2.037 2.033 
(0.03) (0.05) (0.240) (0.48) 

Daily 1.511 1.506 1.512 2.482 2.485 2.482 
(13.96) (37.93) 

Weekly 3.008 2.997 2.963 3.006 3.012 2.989 
(2.80) (9.47) 



The fixed options parameters are S = 100, T = 0.5, K = 100, cr = 0.20, 
r = 0.10. For the down-and-out, B = 99.9, while for the double knock-out, 
Bi = 95, -B„ = 125. C-V denotes results obtained in [8] while Zvan denotes results 
obtained results obtained in 30, ■ Numbers within parentheses are CPU times 



options to a Wiener-Hopf integral equation, and given a formal inverse z-transform 
solution to this equation in terms of a special function plus infinite sums of simple 
functions. This however does not give rise to exact option values, as the solutions 
need to be evaluated numerically. In addition to the Wiener-Hopf (wh) method 
of [12] the other numerical methods which are compared with the HOBIS method 
in Table [5] are the Markov Chain method (MCh) of [TD], the trinomial tree (tt) 
of [5], the Simpson recursive quadrature method (sq) of [13], and the Monte Carlo 
simulation method (mc) of [3] with antithetic variables and 10® simulation runs. 
Table [5] shows that the hobis values are much closer to the WH values of [T^ which 
is the only numerical method derived from an analytical solution. The table also 
shows that apart from the tt values in column 6, there is in general a good level of 
agreement between all the numerical methods. It should be noted that values from 
the TT method that seems to perform relatively poorly in the table are exactly the 
same values reported for the same discretely monitored barrier option in [5] , some of 
which values appear in Table [5] above. This shows that for discrete monitoring, the 
discrepancies between the hobis values and other values listed in Table [5] cannot 
be perceived as a lack of performance of the hobis. In fact Table [6] simply confirms 
a wide spread perception that tree methods are less efficient compared to the most 
common numerical methods for option pricing. The fixed option parameters in 
TablelUare S" = 100, T = 0.5, X = 100, ct = 0.20, r = 0.10, and iV and B 
denote as usual the number of monitoring dates and the barrier value, respectively. 

Analytical pricing formulas for double barrier options with rebate payment has 
recently been obtained in [25] based on the inversion of the Laplace transform of the 
probability density functions by contour integration. Sidenius obtained a similar 
solution in [25] using an approach based on path counting. These two results are 
certainly among the first ones incorporating a rebate payment for continuous barrier 
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Table 6. Comparison of the HOBIS results for a discretely monitored 
down-and-out call. The parameter set is 5 = 100, T = 0.5, K = 100, a = 
0.20, r = 0.10, and N and B denote as usual the number of monitoring dates 
and the barrier value, respectively. 



N 


B 


HOBIS 


WH-IR 


MCh 


TT 


SQ 


MC 


25 


95 


6.63176 


6.63156 


6.6307 


6.6181 


6.6317 


6.63204 


25 


99.5 


3.35542 


3.35558 


3.3552 


3.3122 


3.3564 


3.35584 


25 


99.9 


3.00848 


3.00887 


3.0095 


2.9626 


3.0098 


3.00918 


125 
125 
125 


95 

99.5 

99.9 


6.16797 
1.96143 
1.51098 


6.16864 
1.96130 
1.51068 


6.1678 
1.9617 
1.5138 


6.1692 
1.9624 
1.5116 


6.1687 
1.9628 
1.5123 


6.16879 
1.96142 
1.5105 



options, and actual barrier option values computed with these formulas hardly 
exists, partly because any computation using these formulas is still quite computer 
intensive. The scarcity of numerical option values in the scientific literature for 
discretely applied barrier options is essentially due to the lack of exact solutions for 
such options. We have thus listed in Table[7]some numerical values for a number of 
knock-out option values that incorporate dividend {q) and rebate {Rb) payments. 
The fixed options parameters are 5* = 100, T = 0.5, K = 100, a = 0.20, r = 0.10. 
The various knock-out option types considered are the down-and-out, the up-and- 
out, and the double knock-out barrier options. The barrier application is either 
continuous or discrete in time. Numbers within parentheses in the table are CPU 
times. 

In Table [7] the values of q and Rb are generally chosen to display some of the 
effects of these parameters. For example if we denote by Cbar and Dbar the value 
of a given barrier option under continuous monitoring and discrete monitoring 
respectively, and by Van the value of the corresponding vanilla option with same 
parameters, in the absence of rebates we must have 

Cbar < Dbar < Van (5.2) 

and the generally very slow convergence of a discretely monitored barrier option 
can be verified using the inequalities in (|5.2p . given that the absolute error in the 
computation of Dbar is at most the difference (Van — Cbar), and this yields a 
better approximation of Dbar, the more {Van — Cbar) is small. Equation (|5.2p is 
no longer true when the option pays a rebate, and many such examples appear in 
the table. The table also confirms that no matter what the options' monitoring 
frequency, dividend payments tend to drag down the call option value. 

As far as the computing times are concerned. Table [7] shows that the CPU times 
in seconds for all the three types of knock-out options considered and for the given 
parameters are almost all between and 0.5 seconds under continuous monitoring 
and do not exceed 1.18 seconds in any case. The computing times under continuous 
monitoring are smaller for values very close to the barrier, as already observed in 
Table [H As for discretely monitored options, the computing time to achieve a very 
high accuracy can reach five thousand seconds, partly because it would be time 
consuming to try to find out the minimum amount of CPU time required in such 
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Table 7. Various HOBIS option values for continuously monitored and dis- 
cretely monitored knock-out calls that incorporate rebate and dividend pay- 
ments. The fixed parameters are S = 100, T = 0.5, K = 100, a = 0.20, r = 
0.10. For discrete monitoring the number of monitoring periods is 25 and the 
mesh size used is 4000 X 6000. Numbers within parentheses are CPU times. 



Bi Bu q Rb Continuous Discrete 







Down - and - Out 




90 




0.05 


3 


7.491 (0.42) 


7.493 


90 




0.05 


1.125 


6.719 (0.55) 


6.841 


99.9 




0.05 


3 


3.106 (0.04) 


4.924 


99.9 










0.165 (0.03) 


3.009 








Up - and 


-Out 






110 


0.20 





0.225 (1.18) 


0.344 




110 


0.02 





0.299 (1.15) 


0.470 




100.1 


0.0 


0.01 


0.009 (0.02) 


0.009 




100.1 


0.0 


3 


2.987 (0.04) 


2.735 






Double Knock - Out 




95 


125 








2.033 (0.240) 


3.008 


95 


125 


0.04 


6.66 


7.057 (0.30) 


7.256 


75 


185 


0.045 





6.863 (0.30) 


6.864 


80 


120 


0.04 





2.196 (0.4) 


2.654 



case of very slow convergence to achieve a given level of accuracy. We have however 
noted that the computation of the discretely monitored double knock-out option 
with q = 0.045 in the table gives the indicated value of 6.864 in only 4.20 seconds 
which turns out to be the limiting value, and this is truly a record CPU time for 
the computation of a discretely applied barrier option. 

All of our nonuniform grid implementations as described in Section |4] did not 
yield any faster convergence or any better accuracy in the transformed {x,t) - 
coordinates . As in Section [H denote by xj the jth grid point in the x - direction, 
and by hj — Xjj^-i — Xj the jth grid increment, and similarly let Sj = Ke^^ and h* 
be the corresponding grid point and grid increment in the original S'-coordinates. 
Then we readily see that 

showing that a uniform grid in the x-coordinates (corresponding to 
hj = h = constant) always yields a nonuniform grid in the original S'-coordinates. 
Moreover, such a nonuniform grid is typically sufficiently refined, since h is assumed 
to be small and thus e'* — l is close to zero. For instance, if ft. = 1/10^ and Sj = 50, 
then h* = 0.00500025. This shows that the implementation of effective nonuniform 
grids is not straightforward under the combination of the coordinates transforma- 
tion and the high-order accurate FDS. For all of our attempts of nonuniform grids 
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in the x- or r-coordinates, we achieved the best results only when the non unifor- 
mity parameter qj = hj/hj-i was equal to 1, which corresponds to the uniform 
grids discussed in this section. 

6. Conclusions 

For the valuation of various barrier options, wc have used in this paper a finite 
difference scheme which is fourth-order accurate in the space variable and second- 
order accurate in the time variable. This scheme has been shown to be much faster 
and accurate then all of the most common schemes, namely the Crank- Nicolson, the 
explicit, and the fully implicit schemes. We've also given an optimal determination 
of the boundary conditions, which is particularly relevant for single-barrier options, 
but also for discretely monitored double knock-out options. In this way, with a 
reasonable accuracy we have been able to value continuously monitored barrier 
options in a fraction of seconds and with a mere 20 x 20 mesh. 

The efficiency of the combined high-order scheme and optimal determination of 
boundary conditions was again demonstrated in the difficult problem of valuing 
discretely monitored barrier options. Indeed, a comparison of the option values 
obtained using the resulting scheme with those obtained using five other numerical 
methods shows not only a good agreement, but also that our values are much closer 
to the sole method based on an analytical solution for discretely sampled barrier 
options. Our optimal determination of boundary conditions naturally has several 
applications, and could be used amongst others to classify double barrier options 
into simpler types of options. 
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